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ABSTRACT 

We study the fourth-order moment of the cosmic shear field using the dark matter halo ap- 
proach to describe the nonlinear gravitational evolution of structure in the universe. Since the 
third-order moment of the shear field vanishes because of symmetry, non-Gaussian signatures 
in its one-point statistics emerge at the fourth-order level. We argue that the shear kurtosis 
parameter SL 4 = (jf)c/ {l?} 3 ma y be more directly applicable to realistic data than the 
well-studied higher-order statistics of the convergence field, since obtaining the convergence 
requires a non-local reconstruction from the measured shear field. 

We compare our halo model predictions for the variance, skewness and kurtosis of lens- 
ing fields with ray-tracing simulations of cold dark matter models and find good agreement. 
The shear kurtosis calculation is made tractable by developing approximations for fast and 
accurate evaluations of the 8-dimensional integrals necessary to obtain the shear kurtosis. We 
show that on small angular scales, 9 ^ 5', more than half of the shear kurtosis arises from cor- 
relations within massive dark matter halos with M £ 1O 14 M . The shear kurtosis is sensitive 
to the matter density parameter of the universe, £l m o, and has relatively weak dependences 
on other parameters. Therefore, a detection of the shear kurtosis can be used to break de- 
generacies in determining fi m o and the power spectrum amplitude erg, so far provided from 
measurements of the two-point shear statistics. The approximations we develop for the third- 
and fourth-order moments allow for accurate halo model predictions for the 3-dimensional 
mass distribution as well. We demonstrate their accuracy in the small scale regime, below 2 
Mpc, where analytical approaches used in the literature so far cease to be accurate. 

Key words: cosmology: theory — gravitational lensing — large-scale structure of universe 



1 INTRODUCTION 

Weak gravitational lensing caused by the large-scale structure of the universe has been established as a useful probe of cosmological param- 
eters and offers the possibility of directly measuring the dark matter power spectrum (see Mellier 1999 and Bartelmann & Schneider 2001 
for reviews). Several independen t groups have reported significant detections of lensing by large-scale structur e on distant galaxy images 
(cosntic^hea^jromjhegcom Wilson & Luppino 

200 0; |Maoli et al. 200 It Van Waerbeke et al. 2001a|; Hoekstra et al. 200$ |Bacon et al. 2002|) and from space ( Rhodes, Refregier & Groth 
2001; Haemmerle et al. 2002| ; Refregier, Rhodes & Groth 2002| ). These groups measured the two-point correlation function of the cosmic 
shear field or the variance of the filtered shear field and set constrains on cosmological parameters, in particular some combination of the 
overall amplitude of matter power spectrum (as) and the matter density parameter of the universe (f2 m p), as shown in earlier theoretical work 



dBlandford et al. 199 1| ; |Miralda-Escude 199 1| ; Raiser 1992 
Kaiser 1998|) 



Villumsen 1996 



Bernardeau, Van Waerbeke & Mellier 1997; Jain & Seljak 1997 



It has been shown that the non-Gaussian signature in the weak lensing field induced by nonlinear gravitational clustering can be used to 
break degeneracies in the determination of erg and Q mu (Bernardeau, Van Waerbeke & Mellier 1997; Jain & Seljak 1997). This possibility 
is attractive, since it can determine Q m o via weak lensing measurements without invoking any other methods such as the cosmic microwave 
background (CMB) and galaxy redshift surveys. This also indicates that the dark energy component of the universe can be constrained 
by combining lensing measurements with the evidence for a flat universe revealed by the measured CMB angular power spectrum (e.g., 
Netterfield et al. 2001). Most theoretical work so far has focused on the non-Gaussian signatures described by the higher-order moments of 
the filtered convergence field (Bernardeau et al. 1997; Jain & Seljak 1997; Hui 1999; Van Waerbeke, Bernardeau & Me llier 1997; Jain, Seljak 
& White 2000; Van Waerbeke et al. 2001b; Munshi & Jain 2001) or the skewness parameter in the aperture mass map (Schneider et al. 1998 
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Bartelmann & Schneider 2001 ). It was recently als o proposed that the genus curve or Minkowski functionals of the convergence field could 
be an efficient measure of the non-Gaussian signal ( Sato et al. 2001 ; Matsubara & Jain 2001 ; Taruya et al. 2002 ). Unfortunately these methods 
turn out to have limitations in application to realistic data. Because realistic data has a non-trivial survey geometry with many masked areas 
due to light scattering, bright stars and so on, it is very challenging to reconstruct the convergence from the measured shear field. On the 
other hand, although the aperture mass method has the advantage of being directly obtained from the shear map, it is likely to suffer from low 
signal-to-noise ratio for the skewness measurement, because this method uses a compensated filter (Schneider 1996; |Schneider et al. 1998 ) 
and thus leads to the loss of the non-Gaussian signal, especially on angular scales smaller than a few arcminutes where the signal is large 
(see Van Waerbeke et al. 2001a for detailed comparisons between various two-point statistical measures of the shear field for actual data). 

Very recently, Bernardeau, Van Waerbeke & Mellier (2002a; hereafter BvWM) proposed that some specific patterns in the three-point 
function of the shear field can be used to extract the non-Gaussian signal. Bernardeau, Mellier & Van Waerbeke (2002b) then reported a 
detection of this signal from actual data on 2 — 4 arcminutes scales, although the signal-to-noise so far is not enough to put robust constraints 
on Qmo- The method proposed by BvWM appears to be a promising new measure of non-Gaussianity. It is possible that their method 
loses some non-Gaussian information because the vector-like property of their statistic leads to partial cancellations between the signal on 
averaging. Their method also seems to have the limitation that it cannot extract the signal on small scales (8 S ^ 2'), since their three-point 
function decreases for smaller scales and approaches zero at zero separation. 

The purpose of this paper is to develop an alternative statistical method directly applicable to the cosmic shear data. The method we 
propose is the connected part of the fourth-order moment of the filtered shear field, in particular the shear kurtosis parameter defined by 
S-yA = (jt)c/ {it) 3 , since the non-Gaussian signal appears first at the fourth-order level for the one-point statistics. The kurtosis parameter 
collapses the information in the trispectrum into a single less noisy quantity, although it does not retain the full information in the four-point 
statistics. 

Since the shear field on relevant angular scales is affected by the nonlinear regime of structure formation (e.g., see Jain & Seljak 1997), 
we need a model to correctly describe the redshift evolution and statistical properties of gravitational clustering up to the four-point level. 
The perturbation theory well st udied in the literature may not be adequate for this task. On the one hand, it is known that the so-called 'hyper- 
extended perturbation theory' (Scoccimarro & Friemann 1999) can describe the strongly nonlinear clustering regime (see for Hui 1999, van 
Waerbeke et al. 2001b and Munshi & Jain 2001 for applicatio n to weak lensing). H o wever, the model does no t describe the intermediate-scale 
transition between the linear and strongly nonlinear regimes ( Cooray & Hu 2001a ; Scoccimarro et al. 2001 ), which does affect weak lensing 
statistics on a range of scales because of projection effects. We therefore choose to employ the dark matter halo approach, where gravitational 
clustering is described in terms of correlations between and within dark matter halos (see McClelland & Silk 1977; Peebles 1980; Scherrer 
& Bertschinger 1991 for initial applications; for recent developments see e.g. Sheth & Jain 1997; Komatsu & Kitayama 1999; Seljak 2000; 
Ma & Fry 2000; Scoccimarro et al. 2001; Cooray, Hu & Miralda-Escude 2000; Cooray & Hu 2001a,b; and Cooray & Sheth 2002 for a recent 
review). There are several reasons we use the halo model. First, the halo model is formally complete and simple enough that higher-order 
statistics of the weak lensing fields can be analytically calculated. Second, the results can be interpreted in terms of halo properties, which 
is convenient for comparison with other observations such as X-ray and optical surveys of clusters of galaxies. Finally, the model appears 
remarkably s u ccessful in that, e ven though it relies on ra ther simplified assumptions, it has reproduced results from numerical simulations 
( jScljak 2000| ; |Ma & Fry 2000|; [Scoccimarro et al. 200l[) and also allowe d for interpretations of observational results of galaxy clustering 



( Seljak 200C ; Scoccimarro et al. 2001 ; Guzik & Seljak 2002; Seljak 2002). 



Once the three ingredients of the halo model (halo profile, mass function and halo bias) are specified, it is straightforward to develop 
the formalism to calculate the shear kurtosis. Cooray & Hu (2001a) have investigated the bispectrum of the convergence field using the halo 
model and find the convergence skewness is mainly due to rare and massive halos on relevant scales of 8 S ^ 10', which is referred to as the 
1-halo term in this paper. We will also find that the shear kurtosis arises mainly from the 1-halo term on the relevant scales. However, since 
the direct application of the halo model requires an 8-dimensional integration to obtain the 1-halo term, we develop an approximation that 
significantly reduces the computational time and gives the shear kurtosis with 10% accuracy at most on angular scales of interest. Our model 
predictions will be compared in detail with ray-tracing simulation results for all the statistical measures we investigate: the convergence 
or shear variance, the convergence skewness and the kurtosis parameters of the convergence and shear fields. This comparison addresses 
the broader issue of whether the halo model can accurately describe statistical properties of weak lensing fields for higher-order moments 
beyond the two-point statistics well studied in the literature. We will pay special attention to the dependences of the shear kurtosis on the 
cosmological parameters, f2 m o and as, for flat CDM (cold dark matter) models. 

The outline of this paper is as follows. In §^ we present the dark matter halo model used in this paper and then write down the expressions 
for the power spectrum, bispectrum and trispectrum for the underlying density field. In we investigate the validity of the halo model for 
weak lensing statistics by comparing the predictions with the ray-tracing simulation results for the variance and skewness of the filtered 
convergence field. In we develop an approximation for calculating the convergence kurtosis and extend it to the shear kurtosis calculation 
in The dependence of the shear kurtosis on cosmological parameters is presented in Finally, §j7]is devoted to a summary and discussion. 
In the following, without explicit mention we will often consider two CDM models: one is the SCDM model with f2 m o = 1, h = 0.5 and 
erg = 0.6 and the other is the ACDM model with fi m o = 0.3, Q\o = 0.7, h — 0.5 and as = 0.9, respectively. Here, f2 m o and Q\q are 
the present-day density parameters of matter and cosmological constant, h is the Hubble parameter, and as i s the rms mass fluctuation s of a 
sphere of 8/i _1 Mpc radius. The choice of as for each model is motivated by the cluster abundance analysis (Eke, Cole & Frenk 1996). 



2 DARK MATTER HALO APPROACH 
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2.1 Ingredients 

In the dark matter halo approach the underlying density field can be described in terms of correlations between and within dark matter halos, 
which are taken to be locally biased tracers of density perturbations in the linear regime. The method is therefore based on three essential 
ingredients well studied in the literature: the mass function of dark matter halos, the halo biasing function, and the halo density profile. 

For the halo mass function, we adopt an analytic al fitting model proposed by Sheth & Tormen (1999), which is more accurate on cluster 



mass scales than the original Press-Schechter model (Press & Schechter 1974). The number density of halos with mass in the range between 
M and M + dM is given by 

-JT7 dM = S/M^ 

dM M 

= gA[l + M-V^exp(-f)^, (1) 
where v is the peak height defined by 

5 c (z) 



D(z)a(M) 



(2) 



po is the mean cosmic mass density today (we use comoving coordinates throughout) and the numerical coefficients a and p are empirically 
fitted from N-body simulations as a = 0.707 and p = 0.3. The coefficient A is set by the normalization condition J dvj(y) — 1, leading 

to A ~ 0.129. Here <r(M) is the present-day rms fluctuations in the matter density top-hat smoothed over a scale Rm = (3M/47rpo) 1/ ' 3 , 
D(z) is the growing factor (e.g., see Peebles 1980), and S c (z) is the threshold overdensity for spherical collapse model (see Nakamura & 
Suto 1997 and Henry 2000 for useful fitting functions). It should be noted that the peak height v is given as a function of M at any redshift. 

Mo & White (1996) developed a useful formula to describe the bias relation between the dark matter halo distribution a nd the underlying 
dens ity field. This idea has been improved by several authors using N-body numerical simulations (|Mo. ling & White 1997j; Sheth & Lemson 
1999; |Sheth & Tormen 199^ ; we will use the fitting formula of Sheth & Tormen (1999) for consistency with the mass function (|l|): 

, / % , oxi — 1 2p 

b{u) = 1 + — ■=— + YTTT7 w 0) 

where we have assumed scale-independent bias and neglected the higher order bias functions (62 , 63 , • • •) that have a negligible effect on our 
final results. 

The density profile of dark matter halos is defined to be an average over all halos with a given mass M and does not necessarily assume 
all halos have the same profile and spherical symmetry as stressed by Seljak (2000). It is not evident that this argument should be valid for the 
higher-order moments of the density field or the weak lensing field. However, the agreement between our model predictions and numerical 
simulations indicates that there is no strong violation of the assumption. Throughout this paper we assume the NFW model for the averaged 
halo profile (Navarro, Frenk & White 1996, 1997; hereafter NFW): 

1 ^ / c3M 1 

p( r \ — 1 . ) (4) 

47rr^ ir cr/r v ir(l + cr/r v i v ) 2 ' 

where / = l/[ln(l + c) — c/(l + c)] and r v i r is the virial radius of the halo. The virial radius can be expressed in terms of the halo mass 
M and redshift z based on the spherical collapse model: M = (4nr^ iT /3)poA(z), where A(z) is the overdensity of collapse given as a 
function of redshift (e.g., see Nakamura & Suto 1997 and Henry 2000 for a useful fitting formula). We have A w 340 for the ACDM model. 
It is worth noting that some studies based o n N-body simulations with higher resolution than in NFW have suggested a steeper slope for th e 



inner profile with p oc r -1 ' 5 at r < r vil /c (Fukushige & Makino 1997; Moore et al. 1998; ling & Suto 200C; Fukushige & Makino 2001) 
whereas the predictions for the outer parts of halos are in agreement with NFW: p oc at r £ r v i r /c. Lensing statistics on angular scales 
of interest are affected more strongly by the outer part of the density profile. Further the outer profile is scaled by the concentration parameter 
c for a given virial radius, so we simply assume the NFW profile and pay close attention to the appropriate choice of c as discussed below. 

To give the halo profile in terms of M and z, we further need to express the concentration parameter c in terms of M and 2; however, 
this still remains somewhat uncertain. The concentration c is theoretically expected to be a weak function of halo mass as given by c = 
cq{M / M*)? , where the normalization is Co ~ O(10) at the present-day nonlinear mass scale M* defined by 5 c (z = 0)/cx(M*) — 1 and 
the slope is (3 ~ — O(10 _1 ). We employ the form motivated by Seljak (2000): 

c(M, ,) = 10(1 " (5) 



where we have assumed the redshift dependence (1 + z)^ 1 as supported by numerical simulations (Bullock et al. 2001). There are several 
reasons we adopt the form (^) for the unknown concentration parameter. As for the slope (3, we assume f3 — —0.2 which is steeper than 
(3 — —0.13 originally proposed by NFW and Bullock et al. (2001). This is motivated by the fact that for the NFW profile (Q) the halo model 
with j3 — —0.2 can better reproduce the well-studied nonlinear matter power spectrum than the model with f3 = —0.13 as shown in Seljak 
(2000; also see Cooray et al. 2000). As will be shown below, our model can also reproduce the simulation results for the higher-order statistics 
of weak lensing fields on relevant angular scales. In this sense, for the purpose of using the halo model to describe the nonlinear gravitational 
evolution, it seems to be appropriate to choose the halo model parameters so that the model can reproduce the matter power spectrum as the 
first step. The choice of the normalization of Co = 10 at M* is supported by N-body simulations (Bullock et al. 2001) and also validated by 
the fact that the form (m is consistent with recent observational results of c ~ O(10) on galactic scales of M ~ 1O 12 M0 obtained from 
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analyses of galaxy rotation curves (limenez, Verde & Oh 2002) and galaxy-galaxy lensing ( |5uzik & Seljak 2002 ; Seljak 2002] ). We will 
discuss in more detail how possible variations in the concentration parameter affect the final results of the shear kurtosis. 
The normalized Fourier transform of the NFW profile (W) is given by 



y(k,M; 



1 

M 



4nr 2 dr p(r 



sin kr 
kr ' 



(6) 



where y(k) has the asymptotic behavior y(k) w 1 and y(k) oc k 2 for kr vlr /c <C 1 and kr vlI /c 3> 1, respectively. 



2.2 The power spectrum, bispectrum and trispectrum 

The power spectrum P(k), bispectrum B(k\, fc 2 , fc 3 ) and trispectrum T(ki, fc 2 , fc 3 , fc 4 ) of the dark matter density fluctuation are defined 
by 

(5(ki)S(ka)) = (27r) 3 P(fci)fo(fci 2 ), 
(5(k 1 )6(k 2 )S(k 3 )) = {2-KfB{k 1 ,k 2 ,k 3 )8 D {k 12Z ) 

{5{k{)...5{k 4 )) c = (27r) 3 r(fe 1 ,fe 2 ,fe3,fe 4 )5i 3 (fe 12 34), (7) 

where ki...j = k% + . . . + kj, fo(fe) is the delta function, and (. . .} denotes the ensemble average. The subscript c denotes the connected 
part; the trispectrum is identically zero for a Gaussian field. 

In the picture of the halo approach, the power spectrum can be expressed as the sum of correlations within a single halo (denoted the lh 
term) and between different halos (the 2h term); 

P(k) = P lh (k) + P 2h (k), (8) 
with 

P lh (k) = I°(k,k), 

P 2h (k) = [/!(£)] 2 P L (fc)- (9) 
In the above equations have used the notation of Cooray & Hu (2001a): 

j£(fe,...,fc M ) = ydM^0Q bp{M)y(k 1 ,M)...y{k tl ,M). (10) 

Note that we set bo = 1, &i = 6 given by equation <|3j> and b t = for i > 2. The quantity P L (k) denotes the linear power spectrum, 
and its redshift evolution is given by P L (k,z) = D 2 (z)P L (k,z = 0), although we will often omit z in the argument for simplicity. The 
requirement that on large scales (k — > and y ~ 1) the 2-halo contribution to the power spectrum reduce to the linear power spectrum 
imposes the condition Jduf(v)b(u) — 1, which is automatically satisfied by equations ([!]) and ^ within a few percents. 
Similarly, the bispectrum can be expressed as sum of the 1-halo, 2-halo and 3-halo contributions: 

B — B lh + B 2h + B 3h , (11) 

with 

B lh = Il(ki,k 2 ,k 3 ), 

B 2h = P i (fci)/ 2 1 (fc 2 ,fc 3 )/i 1 (fci)+2perm., 

B 3h = 73 pt (fc 1 ,fc 2 ,fc 3 )/ 1 1 (fci)J 1 1 (/c 2 )J 1 1 (fe 3 ), (12) 

where B pt denotes the bispectrum calculated by perturbation theory and the explicit expression is given in Appendix ^. 
Finally, the trispectrum arises from four contributions involving one to four halos (Cooray & Hu 2001bl: 

T = T lh + (T™ + T 2 2 2 h ) + T 3h + T 4h , (13) 

with 

T lh = lJ(fa,fc»,*B,fc4). (14) 

T31 = P L {k 1 )I 1 i {k 2 ,k 3 ,k 4 )ll{k 1 ) +3 perm., (15) 

TH = P i (fci 2 )/ 2 1 (fc 1 ,fc 2 )/ 2 1 (fc 3 ,fc 4 ) + 2perm., (16) 

T 3h = B pt (fe 1 ,fe 2 ,fe 3 )J 2 1 (fc 3 ,fc 4 )/i 1 (fci)/i 1 (fc 2 ) +5 perm., (17) 
t ah = T Pt( fel) ... jfe4 )/ 1 i( fcl ).../i( A;4 ) ) (1 8) 

where T pt denotes the trispectrum given by perturbation theory (see Appendix Note that the 2-halo term is further divided into two 
contributions, T 31 h and T22 > which represent taking three or two points in the first halo and then one or two in the second halo. 
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3 VALIDITY OF THE HALO MODEL FOR WEAK LENSING STATISTICS 

In this section, we investigate the validity of the halo model to compute weak lensing statistics by comparing our model predictions to 
ray-tracing simulations for the variance and skewness of the filtered convergence field. 



3.1 Weak lensing convergence and shear fields 



The weak lensing convergence is expressed as a weighted projection of the density fluctuation field between source galaxy and the observer 
(e.g., see Mellier 1999 and Bartelmann & Schneider 2001): 

(19) 



k(0) = JdxW(x,Xs)5[x,dA(x)0], 
where x is the comoving distance and the function W is the lensing weight function defined by 

W(x, Xs ) = ~V m0 H a ^ 



X) 



(20) 



Here Hq is the Hubble constant (Ho — lOO/i km s 1 Mpc *) and the function (1a(x) is the comoving angular diameter distance. Note that 
throughout we assume all source galaxies are at a single red shift z s for simplicity. The key simplification used in equation ( fl9| ) is the Born 



approximation (Blandford et al. 1991 



Vliralda-Escude 1991 



Kaiser 1992), where the convergence field is computed along the unperturbed 



path. Jain et al. (2000; hereafter JSW) found that it is an excellent approximation for the two-point statistics. Based on this result, we will 
assume that the Born approximation also holds for the higher-order statistics we are interested in. 

A direct observable of weak lensing is the distortion effect on source galaxy images characterized by the two components of the shear 
field, 71 and 72, which correspond to elongations or compressions along or at 45° to a;-axis, respectively. In Fourier space, the shear fields 
71 and 72 are simply related to the convergence field via the relation 



71 (I) - k(l) cos(2<pi), 72 (Z) = «(Z) sin^i), 



(21) 



where Z = /(cos ipi, sin ipj), quantities with tilde symbol den ote their Fourier components and we have employed the flat-sky approximation 



(Blandford et al. 1991 



Miralda-Escude 1991 



Kaiser 19921. Equation ( plj ) shows that ji has a vector-like property. More specifically, for 
example, each shear component could be either positive or negative even around a dark matter halo on the sky, whereas the convergence field 
is always positive. The statistical symmetry of the shear components around is the reason that all odd moments of the shear field vanish. 
Hence the first non-vanishing non-Gaussian signal appears at the fourth-order level for the one-point statistics. 

In practice spatially filtered lensing fields are used in order to reduce the noise contribution due to the intrinsic ellipticities of source 
galaxies. The filtered shear field can be expressed as 



±L^(l)F(l;8 s )e X p d - 6 . 



Throughout this paper, we employ the top-hat filter function with its Fourier transform given by 



(22) 



(23) 



where J\ (x) is the first-order Bessel function. In the following, we will omit the superscript T for the filtered fields of k and 7 for simplicity. 



3.2 Variance and higher-order moments of the filtered convergence field 

The variance of the filtered convergence field can be expressed as a weighted integral of the dark matter power spectrum: 

fZdZ„/, I 



*l(e s ) = (K 2 (e s ))= d x w 2 (x,xs)d- A 2 (x) h^ p ( k = ^x) F ^ 1 ) 



(24) 



This equation is derived by using the Limber approximation (Limber 1954; also see Kaiser 1992) under the flat sky approximation. It should 
be noted that the angular mode I is related to the three dimensional wavenumber as k = 1/d.A- By using the expression in equation (^) for 
P(k), we can compute the convergence variance based on the halo model. 

Figure [l] plots the convergence variance as a function of the top-hat smoothing scale for the SCDM (f2 m o = 1, h = 0.5, as = 0.6) 
and ACDM (f2 m o = 0.3, Q\o = 0.7, h = 0.7, erg = 0.9) models. We fix z s = 1 for the source galaxy redshift. For the linear matter 
power spectrum used in the calculation, we employ a scale invariant spectrum of the primordial fluctuations with the BBKS transfer function 



(Bardeen et al. 1986). The solid and dashed lines show the results of our halo model for the ACDM and SCDM models, respectively, while 
the dot-dashed lines are the predictions of using the Peacock & Dodds (1996; hereafter PD) fitting formula for the nonlinear power spectra. 
The 1-halo and 2-halo contributions are shown by the thin solid lines for ACDM, and one can see that the variance arises mainly from the 
1-halo term on angular scales of 9 S 5', where nonlinear structures play important role to the weak lensing statistics (see, e.g. Jain & Seljak 
1997). The symbols with error bars are the ray-tracing simulation results, where the error in each bin denotes the sample variance for a weak 
lensing survey with an area of 25 degree 2 . The ray-tracing simulation builds on an N-body simulation based on the particle-mesh (PM) code 
and has been kindly made available to us by T. Hamana (for details see Hamana & Mellier 2001 ; hereafter HM). 

It is clear from Figure |l| that the halo model predictions are in good agreements with the PD results as well as with the simulation results 
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9 s [arcmin] 



Figure 1. The variance of the convergence field as a function of the top-hat smoothing scale 9 B . The solid and dashed lines show the halo model predictions 
for the SCDM and ACDM models, respectively, with the source redshift z B = 1. The square and triangle symbols are results from ray-tracing simulation, with 
error bars giving the sample variance for a survey area of 25 degree 2 , calculated from the simulation data of Hamana & Mellier (2001). The thin solid lines 
are the 1 -halo and 2-halo contributions for the ACDM model. The dot-dashed lines are the predictions from the Peacock-Dodds fitting formula. 
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Figure 2. Shown is the comparison of the high resolution simulation results (diamond symbol) of cr 2 (Jain et al.2000) with the results (square) of lower 
resolution data (Hamana & Mellier 2001) for SCDM. The latter data is mainly used for the comparisons with model predictions in this paper as explained in 
the text. The solid and dot-dashed lines are the predictions of the halo model and PD, respectively, as in Figure |l| Note that the error bars for the JSW data 
correspond to the sample variance for a survey area of 2.8 X 2.8 degree 2 . 



for both the ACDM and SCDM models. We have indeed confirmed that for all cosmological models we consider in this paper the halo model 
can reproduce the PD results for <r 2 (# s ) within ~ 5% accuracy on angular scales of interest. This success at the two-point level is partly 
due to our choice (Q) of the concentration parameter for the NFW profile. However, there are slight differences between the predictions and 
the numerical results on small angular scales 9 S £3'. This is possibly due to the lack of the numerical resolution of the ray-tracing data, 
because the higher resolution simulation used by JSW yields more power on such small scales, which gives a better match to the theoretical 
predictions, as explicitly shown in Figure |^ (see also discussions in Taruya et al. 2002 for the resolution of the HM data). We prefer to use 
the HM data for comparison to our model predictions because we can use 40 realizations of simulation data with 25 degree 2 for each CDM 
model in order to correctly estimate the sample variance. Having an adequate number of the realizations is crucial to study the higher-order 
moments especially on large angular scales, 9 S ^ 5', since the higher moments are more sensitive to sample variance. 

In analogy with the second moment, the third-order moment of the filtered convergence field can be expressed in terms of the bispectrum 
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Figure 3. The skewness parameter as a function of 6 S as in Figure [lj The upper panel shows a comparison of the halo model predictions with the simulation 
results, while the lower panel shows contributions from the 1-halo, 2-halo and 3-halo terms for the ACDM model. For comparison, the dotted line in the lower 
panel shows the result predicted by the second-order perturbation theory. 



<ac 3 (0 s )> = d X W\ X ,Xs)d- A \x) 



d 2 h d 2 l 



- 2 B (fei, fe 2 , -few) F(h)F{l 2 )F{l 12 ) 



(2tt) 2 (2tt) 2J 

where ki — li/dA( X ) an d B f° r the halo model is given by equation dill). We explicitly write down the 1-halo contribution to (k 3 }: 



<K 3 (0 8 )r= /d X W/ 3 (x,X s )dI 4 (x) dM 



dn (MY f d 2 h d 2 l 2 



dM \p 



(2tt) 2 (2tt) ; 



■y(h, M)y(l 2 , M)y(h 2 , M)F(h)F(l 2 )F(l 12 ). 



(25) 



(26) 



Although the Fourier transform of the NFW profile, y(k, M), is given as a function of the three-dimensional wavenumber k, we will often 
use I for the argument of y(k, M) according to the relation of k = l/dA{x) f° r simplicity. To obtain (« 3 (0 S )) , we need to perform a 
5-dimensional numerical integration, since we can eliminate one angular integration using statistical symmetry. The convergence skewness 
parameter is defined by 



(27) 



This form is motivated by the fact that in perturbation theory both the numerator and denominator in equation ( p7[ ) scale as ~ 5f, where 8i 
is the linear solution for the density fluctuation field. Hence the skewness becomes almost independent of the power spectrum normalization 
as, giving roughly a dependence as S K ,3 oc fi~Q through the dependences of the angular distances and the growth rate of the fluctuations 
(Bernardeau et al. 1997). In the results shown below, we use the halo model self-consistently to compute o K (6 s ) in the denominator of S K ,3. 
Since our halo model can reproduce the PD results for a 2 within 5% accuracy, this does not significantly affect our results for the skewness 
or kurtosis parameters. 

Figure ^ plots the convergence skewness parameter as a function of the smoothing scale as in Figure |l|. It is clear from the upper panel 
that the halo model prediction agrees well with the simulation result for SCDM over all scales. For the ACDM model our model slightly 
overestimates the simulation result on small scales. Among the possible reasons for this discrepancy, one is that the HM simulation result 
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Figure 4. The left panel shows the dependences of the integrand on redshift for the variance (solid line) and the third- (dashed) and forth-order (dotted) 
moments of the convergence field. The source redshift and the filtering scale are z s = 1 and 8 B = 1 '. The bold and thin lines are results for the SCDM and 
ACDM models, respectively, and each curve is normalized by the SCDM value at the peak redshift. The figure shows that compared to the SCDM model, the 
amplitude of each integrand decreases for ACDM and the peak redshift shifts to higher z. In the right panel, we show the comoving transverse distance at the 
peak redshift of these integrands against the angular smoothing scale. 



may underestimate the true value of S K ,3 due to a lack of numerical resolution as explained in Figure g[ As shown in the lower right panel 
of Figure 18 in JSW, the high-resolution N-body simulations yields S K) 3 ~ 140 for ACDM on small angular scales, which gives a better 
match to our halo model prediction. However, the precise value of S Kt 3 for the ACDM model in numerical simulations is still perhaps an 
open issue. Independent ray-tracing simulation performed by White & Hu (2000; hereafter WH) indicate S Kt 3 ~ 110 around 8 S = 4'. We 
found that an important difference in the WH simulations is the values of cosmological parameters, since they use as = 1.2 and V = 0.2 
for the ACDM model, whereas JSW and HM used as = 0.9 and T = 0.21. For the cosmological models used in WH, our halo model 
predicts a ~ 15% decrease of S K ,-j at 8 B < 5' compared with Figure ^ and the resulting S Ky s is then marginally consistent with the result 
shown in Figure 9 of WH on the angular scales we have considered. This is probably due to the increase of as from 0.9 to 1.2, which 
affects the skewness in a complex way since nonlinear contributions are significant in both its numerator and denominator. Thus on these 
small scales the expectation from perturbation theory that S K .;i is independent of erg is not exactly valid. The halo model predicts that the 
skewness and kurtosis of lensing fields slightly decrease with increasing as as shown in Figure [Tv| for the shear kurtosis. Another difference 
between the N-body simulation codes used in HM, WH and JSW is that the JSW data is based on the adaptive particle-particle/particle-mesh 
(AP 3 M) N-body simulations (see Jenkins et al. 1998 in more detail), while the HM and WH data are based on the particle-mesh (PM) 
simulations. The AP 3 M method is expected to achieve higher resolution than PM meth od for similar m esh resolution. We were able to use 
a new high-resolution simulation performed with an AP 3 M code using 512 3 particles JHamana 2002 ) to compute the skewness. The new 
data gives S K ,3 w 137 for 8 B = 1', which agrees with the halo model result in Figure pi but we also find S K ,3 w 125 for 9 S = 4', a value 
higher than the simulation result shown in Figure ^, but still lower than the halo model prediction at the \-a level. Finally, we note that S K ,z 
has a stronger dependence on Q m o with decreasing f2 m o; for example, slight decrease of Afl m o = —0.05 leads to a significant change of 
AS«,3 » 21 at 6 S = 1' for flat CDM models with the ACDM model taken to be the fiducial model, while the skewness for the SCDM 
model is almost unchanged with AS K> 3 ~ 1. These scalings are roughly consistent with the expectation scaling S K ,s oc fl~ . Hence, the 
discrepancy between our model and the simulation result for the ACDM model corresponds to a relatively small change in fi m o. 

In the lower panel of Figure [| the 1-halo, 2-halo and 3-halo contributions to S Kt 3 are separately plotted for ACDM . Notice that, 
for example, '1-halo' here means the convergence third-order moment in S K ,3 includes only the 1-halo contribution, but the convergence 
variance used includes the total contributions from the 1-halo and 2-halo terms. It is apparent that the 1-halo term dominates over all scales 
shown and, in particular, contributes ~ 80% of the total at the smallest scale 8 = 1'. This holds for the SCDM model also. The result 
thus hnpTie^MhaiUhe^fij*h^^ et al. 

2001; pooray & Hu 200 1 a| ; pooray & Hu 2001bj ). The dotted line shows the skewness calculated by perturbation theory for the same ACDM 
model, and it significantly underestimates S Kt 3, since the weak lensing field on relevant scales is affected by strongly nonlinear gravitational 
clustering. These features can be more explicitly explained in Figure H. In the left panel, we plot how the integrand functions for the variance 
(solid line) and the third- (dashed) and fourth-order (dotted) moments of the convergence field depend on redshift for z s = 1 and 8 S = 1'. 
Note that the fourth-order moment is computed using the approximation developed below. The bold and thin lines are the results for the 
SCDM and ACDM models, respectively, where each curve is normalized by the SCDM value at the peak redshift. One can readily see that 
the higher-order moments are more sensitive to lower redshift structures and, compared with the result for SCDM, the amplitude of each 
integrand decreases for ACDM and the peak redshift shifts to higher z. The right panel plots the comoving transverse distance at the peak 
redshift of the integrand function as a function of the smoothing scale 8 S , which is defined by A = d J 4(z P cak)6 s . Again, the figure clarifies 
that the higher-order moments are more sensitive to structures on smaller scales; for example, by comparing the solid and dotted lines one 
finds that the transverse scales for the fourth-order moment are smaller than those for the variance by factors of 0.8 and 0.4 at 8 S = 1' and 
10', respectively 



© 0000 RAS, MNRAS 000,[I|-|^ 



The Kurtosis of the Cosmic Shear Field 9 



It is worth noting differences between the convergence skewness, the shear three-point correlation function recently proposed in BvWM 
(see also Bernardeau, Mellier & Van Waerbeke 2002), and the shear kurtosis. As shown in Figure ||, S K ,3 has a weak dependence on the 
angular scales as pointed out based on the perturbation theory (Bernardeau et al. 1997), while the shear three-point correlation function has 
a logarithmically decreasing behavior with decreasing the angular scale as shown in Figure 6 in BvWM. It is likely that their shear three- 
point function loses useful non-Gaussian information resulting from cancellations between signals caused by the vector-like property of the 
shear field. An advantage of the shear kurtosis parameter is that it collapses information from the 4-point statistics into a single quantity 
without being affected by such cancellations. However the kurtosis parameter is a higher-order moment, and so it remains to be seen how its 
signal-to-noise properties compare with the three-point function of BvWM. 



4 THE KURTOSIS OF THE CONVERGENCE FIELD 

In this section, we develop a useful approximation for fast and accurate evaluations of the convergence kurtosis parameter. In particular, 
we concentrate on developing approximations for calculating the 1-halo term, {k 4 (6 B ))l h , which provides the dominant contribution to the 
convergence kurtosis on small angular scales. The approximations for the 2-halo and 3-halo terms are presented in Appendix Those 
approximations will be used to develop a method to compute the shear kurtosis in the next section. 



4.1 Definition 

The connected part of the fourth-order moment of the convergence field is given by 

{k 4 (6 b )) c = [d X W 4 ( X , Xs)d- A 6 y ^ ^ ^ r(fcr, fc a , fc 3 , -k 123 )F(h)F(l 2 )F(l 3 )F(l 123 ), (28) 

where ki = li/d,A(x) an ^ me trispectrum is given by equation ( fl3| ) within the framework of the halo model. From the results of the 
convergence skewness, it is expected that the most important contribution to the fourth-order moment is the 1-halo term on angular scales of 
interest. The 1-halo term is given by 

("V.))? - Jd X W\ X )d-/ JdM^ (|)7 (2^ ^ y(^M)y(h,M)y{l 3 ,MUl 123 ,M) 

xF(h)F(l 2 )F(l 3 )F(h 23 ). (29) 

Hence, to obtain (K 4 (9 s ))l h , we have to perform at least a 7 dimensional numerical integration, even after we eliminate one angular inte- 
gration of li using statistical symmetry. Direct integration is not suitable for our final purpose of evaluating the dependence on cosmological 
parameters, which requires lots of computations in parameter space. We therefore explore an approximation for calculating (k 4 ) c with 
adequate accuracy and reasonable computational expense. 

Motivated by perturbation theory, as for the skewness parameter, we consider the convergence kurtosis parameter defined by 

We use self-consistently the halo model to calculate a K (9 s ) in the denominator of S K ,4. It is expected that S K ,4 has a dependence roughly 
given by perturbation theory as S K ,t oc fi~g. 

The approximation for calculating S 1 .^ developed below allows us to simplify the three-dimensional angular integrations of 
d 2 lid 2 l2d 2 l3 in equation (^g|), whereby we can obtain (K 4 (6 B ))l h by a 5-dimensional numerical integration instead of the original 8- 
dimensional one. 



4.2 Quadrilateral configuration dependence 

Equation ( p9~| ) shows that, although the integrand function of (K 4 (9 a ))l h does depend on quadrilateral configuration with four sides Zi, I2, I3 
and l±(= — Z123) in Fourier space, the angular dependences of li appear only via Z123 in y(li23, M) and F(l\2i) as a result of the spherical 
symmetry of the NFW profile^. Because of statistical symmetry, without loss of generality we can express any configurations in terms of 5 
parameters; three side lengths of Zr, I2 and Z3 and two angles $2 and $3, where $2 is angle between h and I2 and $3 angle between /12 and 
l 3 . The side length /123 can then be expressed as 

Zi23 = -2Z12Z3 cos $3, (31) 

with Z12 = (I 2 + I 2 — 2Z1Z2 cos 3>2) 1 ^ 2 - Note that the volume element ( p9~| ) of integration can be rewritten, after performing one of the angular 
integrals, as d 2 lid 2 hd 2 l 3 = {2n)hdhl2dl2hdhd$2d<& 3 . 



t On the other hand, when (k 4 (8 s )) c is calculated in perturbation theory, the angular dep ende nces appear via products of the li vectors in the perturbation 
trispectrum in addition to via Z123 in y(/i23) and -F(/i23), as explicitly shown in equation ( |A2| ). 
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4.3 Approximation for the integration of the top-hat filter function 

First, we consider an approximation for the integration of the top-hat filter function (^) motivated by Appendix A in Bernardeau (1994), 
where the geometrical properties of the integration of products of the three-dimensional top-hat window function are derived. In Appendix 
O we prove the following identity for the integration of products of top-hat kernels: 



d 2 h d 2 l 2 d 2 l 
(2tt) 2 (2tt) 2 (2tt) 2 



-I (h)F(l 2 )F(l 3 )F(l 123 ) = ff[ lj ^F 2 (h). (32) 

i — 1 



The result above cannot be applied exactly to simplify equation (g9J) because of the y(li2 3 ) term. We therefore use the following replacement 
for the filter function F(li2 3 ) in equation ( |2^ ) as an approximation to be tested: 

F(li2 3 ) « F(h)F(h)F(l 3 ). (33) 

The corresponding approximation for the three-dimensional window function is used in Scoccimarro et al. (2001) for the study of the 
skewness and kurtosis parameters of the three-dimensional density field. It is worth noting that this approximation indeed becomes exact if 
y(l, M) — constant in equation (^). Hence, to the extent that the regime y(l, M) ~ 1 at I dA(x) c / r vir provides the main contribution 
to (K 4 (9 B ))l h for a given M and z, it is reasonable to expect that the replacement (133) is a good approximation for realistic density profiles. 



4.4 Approximation for the convergence skewness 

Given the approximation (p3|), the next problem we consider is to explore an approximation to describe the configuration dependence of 
y(li23, M) in a way that allows us to evaluate the angular integrations with respect to $2 and $3 in equation (p9|). 

For this purpose, let us begin by considering an approximation for calculating the 1-halo term in the convergence third-order moment, 
(k 3 ) . This is because the accuracy of our approximation for (ft 3 ) 1 ' 1 can be tested by comparing the prediction with the true value obtained 
by the direct integration, and then it can be extended to the calculation of the fourth-order moment. The dependence on the triangle con- 
figuration appears via Z12 in y(li2, M) with I12 = (I 2 — 2l\l 2 cos $ 2 ) 1/2 . We propose a method to expand y(li2) around a fiducial 
triangle configuration with a fixed $2 in analogy with the Taylor expansion of yiliz) with respect to $2, whereby we can analytically per- 
form the angular integrations of ij in equation (p6[). The critical question that arises is: which fiducial configuration is appropriate for the 
expansion? This can be answered by using the halo model analysis of Cooray & Hu (2001a) for the convergence bispectrum, which is part of 
the integrand of (ft 3 ). Figure 7 in their paper explicitly illustrates the configuration dependence of the bispectrum and implies that the main 
contribution to (k 3 ((9 s )) arises from equilateral triangle configurations with l\ = I2 = li2- Hence, it will be reasonable to take a prescription 
that the fiducial configuration contains equilateral triangle configurations when l\ — h. This holds for <E>2 = 7r/3. We thus propose the 
following approximation for calculating {n 3 (6 B )) lh combined with the approximation (j?3]): 

/ \ 3 



<^(0 s )) lft « jd x w - i ( x ,xs)d- A J dM ^[—) J 2 7 ^ y{h)yii2)v{h2)F\h)F\h), 04) 

with l\2 = Qi + I2 — ilia) 1 - Note that I12 = h = I2 when l\ — I2, so that the dimension of integration is reduced from 5 to 4 in 
equation (^rj). Like the Taylor expansion, one can include higher-order corrections arising from the expansion of y(li2, M) at the order of 
0($2 — "Y^Jq We find that the zeroth-order approximation ( p4[ ) works remarkably well as shown below. 

Figure g demonstrates the accuracy of our approximation (^) for the 1-halo term of the convergence skewness by comparing the 
predictions with the direct integration results of equation (^) for the SCDM and ACDM models. The approximation is very accurate, 
as its relative accuracy is better than 5% over all angular scales for both models. For comparison, the dotted and dot-dashed lines show 
the results of using other possible approximations for ACDM, where we used the replacements of y(l\2,M) = y(h, M)y(l%, M) or 
y(h 2 , M) = y(h,M) in equation (Eg), respectively, in addition to the approximation (M) for the filter function. The former approximation 
is motivated on the analogy of equation (|33[), while the latter is indeed used by Scoccimarro et al. (2001) for calculations of the skewness 
and kurtosis parameters of the three-dimensional density field. It is clear that the approximation of y (Z12) = y{h) overestimates the value 
of Sl_ h 3 (see also Cooray & Hu 2001a) and the discrepancy is larger on smaller scales. In more detail, it overestimates S^ 3 by ~ 40% at 
8 B = 1'. The approximation y(lrz, M) = y(h, M)y(l2, M) underestimates the skewness, since y(l, M) < 1, and yields S^ 3 smaller by 
~ 60% than the correct value. Hence, we cannot use these approximations to predict the higher-order moments of weak lensing fields with 
sufficient accuracy for our purpose. 

It is straightforward to apply our approximation to evaluations of the skewness parameter of the three-dimensional density field, = 
{S 3 } j (8 2 ) 2 , which is relevant for surveys of galaxies clustering (e.g., see Scoccimarro et al. 2001). Figure |^ plots the result against the three- 
dimensional smoothing scale R (Mpc) for the ACDM model and 2 = as shown in Figure ^. Note that we have used the three-dimensional 
top-hat filter function, and equation (b3j) can be used as an approximation for the kernel. It is clear that our approximation again works well, 
implying that it will be useful for efficiently exploring parameter space for constraining cosmological parameters from Sg, 3 measurements 
down to very small scales. The issue of the small scale behavior of higher-order moments is somewhat an open question since results from 
numerical simulations are not yet reliable for scales below 1 Mpc. The results shown in this paper for the third and fourth-order moment 



§ In this case, the expansion parameter (<3>2 
guaranteed. 



— 7r/3) could be larger than unity in the range of $2 = 



[0, 2tt], so the convergence of the expansion is no longer 
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Figure 5. Comparison of our approximation for the 1-halo term of convergence skewness parameter with direct integration values for the SCDM and ACDM 
models. The solid lines are the direct integration values, while the dashed lines show the results of approximation (^4|. For comparison, we also show the 
results calculated using other approximations for ACDM: the dot-dashed and dotted lines are computed using the replacements of j/(/i2 , M) m y(li , M) and 
y(Zi2 ) ~ y(h , M)y(l2 , M), respectively, for the integration. In the lower panel, the relative errors of our approximation are shown for the two models. 



are encouraging. To the extent that the current halo model describes clustering accurately, we have tractable analytical means of predicting 
higher-order clustering statistics extending to very small scales. 

4.5 Approximation for the convergence kurtosis 

Based on the above success of our approximation for the skewness parameter, we extend it to develop an approximation for the 1-halo term of 
the convergence fourth-order moment, (k 4 ) \ h . The problem is to consider an efficient expansion of 3/(^123, M) in equation with respect 
to the two angles $2 and $3, such that we can analytically perform the two-dimensional angular integrations of $2 and $3. In the same spirit 
as the approximation (p4|), we choose the fiducial configuration with h — h ~ h ~ ^123, because we believe, in analogy with the skewness, 
that the trispectrum with such configurations produces the main contribution to the fourth-order moment. For the kurtosis, we need to make 
additional choices for the angles $2 and $3; we simply set $2 = n/2 and $3 = 7r/4, which implies a square shaped configuration when 
l\ — h — h- The sketch in Figure ^ illustrates the fiducial 4-point configuration. Applying the above approximation to equation (^) gives 

A/n \\lh ^, I j. .Tir4/. . . . N j-6 /jjix dn M 



<« 4 (0 8 )>r« JdxW *(x,X»)dT JdM— I — J l(M,z;9 s ), (35) 
with 

/3 
Y[ l J^ y ( h ,M)y(h, M)y{h, M)y(l 123 , M)F 2 (h)F 2 (l 2 )F 2 (l 3 ). (36) 

i — 1 

where [i 23 = (^i2 + ^3"V / 2^i2^3) 1 ^ 2 with l\2 = (if + l 2 ) 1 ^ 2 ■ Consequently, to obtain (k 4 )^, we need to perform a 5-dimensional numerical 
integration, which requires much less computational time compared with the original 7-dimensional integration of equation (p^). 

Figure ^ demonstrates the accuracy of our approximation (p5[). The approximate result for T{M, z) is compared with the direct integra- 
tion value, plotted against halo mass M for the ACDM model. The lens redshift is z = 0.4 and smoothing scale 9 S = 1'. Note that for fixed 
2 and M we can directly compute X(M, z) by evaluating the 5-dimensional integral. It is also worth noting that z — 0.4 is chosen because 
it is close to the peak of the lensing weight function W(x, Xs) f° r source redshift z s = 1. The figure clearly shows that for both SCDM 
and ACDM models our approximation can reproduce T(M, z) within 10% accuracy on mass scales of 1O 14 M0 H M ,$ 1O 16 M0, which 
provide the dominant contributions to the kurtosis parameter on relevant angular scales as shown in Figure [l3|. The approximation works 
better for more massive halos. 
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3D Skewness 
ACDM, z = 
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Figure 6. The skewness parameter of the three-dimensional density field, Ss : s = (<5 3 ) /{<5 2 } 2 , as a function of the top-hat smoothing scale, R (Mpc), for the 
ACDM model at z = 0. We demonstrate the performance of our approximation for the 1-halo contribution as shown in Figure ^| The dotted and dot-dashed 
lines show other approximations for the 1-halo term used in the figure, as in Figure |^. The thin lines are the 2-halo, 3-halo and total contributions. 




Figure 7. The sketch of the fiducial 4-point configuration used in the approximation (j^]) for calculation the convergence fourth-order moment. The two angles 
<E>2 and $3 are set to be it/2 and 7r/4, respectively, but the side length parameters ii, li and Z3 are treated as variables. 



To estimate the final accuracy of our approximation to (re 4 )*' 1 , we further need to take into account the lens weighting, W /4 (x, Xs)*^ 6 , 
as well as the weighting of mass function, f(y, z)M 3 , in equation (p5|). Since the lens weighting gives a smooth redshift dependence, we 
here consider the weighing of mass function. The lower panel in Figure ^ plots f(y, z)M 3 at z — 0.4 against M, where each curve is 
normalized to give unity at the peak mass scale. Accounting for the weighting of f(u)M 3 , we find that the accuracy of our approximation is 
about 9% and 8% for SCDM and ACDM models at 6 S — 1' and z — 0.4. Our approximation works better for the larger smoothing scales, 
where more massive halos contribute to (ft 4 ) 4h (see Figure From these results, we are confident that our approximation can predict 
the 1-halo term within ~ 10% accuracy at most on relevant scales, although we should bear in mind that the approximation has a trend to 
underestimate the true value. The figure also shows the results from other possible approximations for ACDM , as in Figure ||, where we have 
used the replacements of y (Z123 ) = y(h)y(h)y(l3) (dotted line) and y{l\2s) = y(lz) (dot-dashed line) (Scoccimarro et al. 2001). These 
approximations overestimate or underestimate I(M, z) by 62% or 70%, respectively, and become worse at more massive mass scales, and 
thus are not accurate enough for our purpose. 

Similarly, we can construct approximations for the 2-halo and 3-halo terms to predict a total power of the convergence kurtosis. The 
explicit forms of the approximations used are presented in Appendix ^. We have confirmed that these approximations are adequately accurate 
(see Scoccimarro et al. 2001 for similar discussions on the skewness and kurtosis of the three-dimensional density field). As explained below, 
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Figure 8. Shown is the accuracy of our approximation ( pq ) for the integrand function X(M, z) for the 1-halo term of the convergence fourth-order moment. 
For the ACDM and SCDM models, the relative errors defined by (X app — X true ) /X true are plotted as the solid and dashed lines, respectively, as a function 
of halo mass M. We here fix z = 0.4 and 8 B = 1' for the lens redshift and the smoothing scale. For comparison, the dotted and dot-dashed lines show the 
results of using other approximations as in Figure ^|(see text for more details). In the lower panel, we plot the mass function weighing, /(i/)M 3 , in the 1-halo 
term for the models, where each curve is normalized to give unity at the peak scale. 



in this paper we ignore the 4-halo contribution which is likely to have a negligible contribution on the angular scales we have considered (see 
Cooray & Hu 2001b for the trispectrum). 

We can now compare our model predictions of the convergence kurtosis parameter with the simulation results. Figure ^ plots the result 
as in Figure Q It is apparent that our halo model predictions are in good agreements with the simulation results as for the skewness case. One 
caveat we should bear in mind is again that the simulation result for ACDM is likely to underestimate S K ,i because of the reasons given for 
Figure^. We have confirmed this by using new high-resolution simulation data provided by Hamana (2002). We obtained S K ,4 ~ 4 x 10 4 
at f? a = 1' for ACDM, which gives a better match to our model prediction. The main cosmological implication of this figure is that there are 
still significant differences between SCDM and ACDM models on small scales of 9 S 2', although the sampling errors corresponding to a 
survey area of 25 degree 2 become larger compared with the skewness case. 

The lower panel of Figure ^ plots the 1-halo, 2-halo and 3-halo contributions for ACDM. It is clear that the 1-halo term gives the 
dominant contribution over the scales considered; the 2-halo is marginally important on larger scales of 8 S ^ 5' and the 3-halo makes only a 
small contribution. More explicitly, these terms provide 82%, 16% and 2% of the contributions to the total kurtosis at 9 S = 1'; 53%, 37% 
and 10% at S = 10', respectively. These results validate our expectation that the 4-halo term is negligible on small angular scales 6 S £ 5', 
and even for larger scales it is likely to have contributions smaller than 10%. The 4-halo term is difficult to evaluate by numerical integration 
because of the oscillatory shape of the perturbation theory trispectrum, resulting from its dependences on the interior angles of the 4-point 
configuration. 



5 APPROXIMATION FOR THE SHEAR KURTOSIS 

Based on the results shown in the preceding section, we develop an approximate method for calculating the shear kurtosis, which is the main 
purpose of this paper. 
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Figure 9. Convergence kurtosis S Ky 4 as a function of the smoothing scale as in Figure |3|. The square and triangle symbols with error bars are the simulation 
results for the ACDM and SCDM models, respectively, while the solid and dashed lines denote our model predictions. For illustration, we slightly shift the 
simulation result for SCDM along the x-axis. The thin solid lines in the lower panel are the 1-halo, 2-halo and 3-halo contributions for the ACDM model. For 
comparison, the dotted and dot-dashed lines are the results of other approximations for the 1-halo term as in Figure pi 



The connected fourth-order moment of the filtered shear field can be expressed in terms of the convergence trispectrum as 

2 

II T^W cos 2<t>kF(h)T K {h, l2,h,h)^fS D (h2Zi), (37) 

with 

T K = [d X W i (x,Xs)dA 6 T(k 1 ,k 2 ,k 3 ,k 4 ) , (38) 



where k = Z/aU(x)> 4>U ls defined by li — /i(cos t/!^ , sin c^) and T K is the convergence trispectrum (see Cooray & Hu 2001b). This 
equation clarifies that the integrand function of (7i) c has configuration dependences via the geometrical factors of cos2<^; i in addition to 
the convergence trispectrum, and we therefore have to consider the 8-dimensional integration. Note that for (72)0 the geometrical factor in 
equation ( p^ ) is Yit=i s ^ n '^'f'k' rjut (7i)c ~ (72)0 from statistical symmetry. In comparing equation ( [37| ) with equation ( j29| ) for (k 4 ) c , the 
difference is only the geometrical factor Ylt-i cos ■ ^ e therefore expect that the following simple relation 

<7i 4 (0 s )>c w / s (k 4 (0 8 )>c (39) 

applies, with a constant factor f g . We can derive an upper limit for f g in the following rough manner. From the integrand function of (7i) c , 
we consider the angular averaged geometrical function as a function of l\, I2 and Z3 defined by 

Q(l-L,h,h)= /jl^ cos 2^ cos 2cj> ll23 , (40) 

i — 1 

where cos 0; 123 = — (h cos 4>h + h cos <f)i 2 + 1$ cos 4>i 3 )/li23- The function Q peaks at l\ — I2 — h and approaches zero for h <C l 2 ,h 
or li 3> I2, 13 and so on, so that the shear fourth-order moment is suppressed compared to that of the convergence. Hence, the upper limit on 
f g should be set by the case h = h = h'- f g < Q(l, I, I) ~ 5.17 x 10 -2 . 
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Figure 10. Shown is the ratio of the shear kurtosis parameter to the convergence kurtosis as a function of the smoothing scale. The symbols with error bars 
are the simulation results for the ACDM and SCDM models as in Figure |l| where the shear kurtosis is obtained by averaging the two kurtosis parameters for 
71 and 72 and the error bar in each bin is properly averaged. The solid line shows the ratio value of 0.263, which is the average value between the results at 
8 B = 1' for the two models. 



Given the relation (p9|), the shear kurtosis can be expressed in terms of the convergence kurtosis as 

s lA m = -MK - s/.i^r^ = *fMW> (41) 

<T~t(V s ) cr K (tf s J 

where <7 7 (6> s ) is therms of the filtered shear field defined by cr T (0 s ) = (7 l 2 (f? s )} 1 ^ 2 and the factor of 8 comes from the relation (7^} = (re 2 }/2. 
Note that the upper limit on f g discussed above corresponds to S 7 ,4 < 0.415 K ,4. 

Unfortunately, it is difficult to analytically derive the constant factor connecting S 7 ,4 and S K ,4, and we therefore rely on ray-tracing 
simulations. Figure |Io|plots the simulation results for the ratio of S-,,4 to S k a as a function of the smoothing scale. Note that we have taken 
the average of the kurtosis values for the two independent shear fields to obtain Sj^; S^ t 4, = [(7i) c / (7i) 3 + (72}c/(7l} 3 ]/2. The figure 
reveals that, despite the fact that the trispectrum has a strong dependence on cosmological models (leading to a difference greater than 300% 
between the kurtosis values for the SCDM and ACDM models - see Figure the ratios are similar for the two models. Further, it appears 
that S 7l 4 is related to S Kz 4 by a constant factor over the angular scales we have considered. The solid line shows the average value of 0.263 
between the results at 8 S — 1' for the two models, and one can see that the curve reasonably explains the simulation results. In the following, 
to predict the shear kurtosis parameter, we simply multiply the factor 0.263 by the convergence kurtosis parameter calculated using the 



approximations developed in §4.5. Thus we use: 

S 7 ,4 » 0.263 S kA . (42) 



6 RESULTS 

In Figure [IT] we compare our model predictions of the shear kurtosis with the simulation results for the SCDM and ACDM models as in 
Figure ^. Note that the halo model prediction is calculated from the sum of the 1-halo, 2-halo and 3-halo contributions to the shear fourth- 
order moment. The result shown in each bin is computed from the average value of the kurtosis parameters for two shear fields, 71 and 72, 
as in Figure [Hj. The figure reveals that our model can well reproduce the simulation results and that there are distinct differences between 
the shear kurtosis values for the SCDM and ACDM models on small scales of 6 B 3'. The range of angular scales with 0.5' 8 S ,$ 3' is 
feasible for making adequate signal-to-noise measurements of top-hat smoothed statistics from lensing survey data (e.g., see Van Waerbeke et 
al. 2001a). As mentioned in the discussion of Figure | 3| the simulation results for ACDM are likely to underestimate S* 7 ,4 at 8 S ,$ 3'. We have 
confirmed that high-resolution simulation ( Hamana 2002 ) does give a better match to our model prediction, but the accurate measurement 
of fourth-order statistics from numerical data needs further investigation. This uncertainty does not seriously undermine our conclusions 
about parameter estimation from the shear kurtosis, because it has a strong dependence on il m o for flat ACDM models. E.g. a small change 
Afi m o = —0.05 leads to a large change of AS 7 ,4 w 3.1 x 10 3 at 6 S = 1' if one chooses the fiducial model with fi m o = 0.3, Q\o = 0.7. On 
the other hand for variations about the fi mU = 1 model, the result is almost unchanged with AS 7 4 w 60. These scalings are approximately 
consistent with the perturbation theory expectation given by S 7 ,4 oc Sl~g (see also Figure [n|and[jj). 

In Figure |l2| we show the comparison of our model prediction with the JSW simulation results for the rCDM model, which has 
fimo = 1-0, h — 0.5, as = 0.6 and the shape parameter of T = 0.21 (see ISW for more details). Note that the cosmological parameters 
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Figure 11. Comparison of our model predictions for the shear kurtosis, 4, with the simulation results as a function of the smoothing scale. Here, for the 
simulation result in each bin we take the average between the kurtosis parameters of two shear fields 71 and 72 . For illustration, we slightly shift the simulation 
result for SCDM along the x-axis. 




Figure 12. Shown is the comparison of our model predictions for the shear kurtosis (upper panel) and the convergence skewness (lower panel) with the results 
of the JSW simulation data for the rCDM model, which has S7 m o = 1.0, h = 0.5, as = 0.6 and V = 0.21. 
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Figure 14. The dependence of the shear kurtosis on the range of lens redshift is plotted. From top to bottom, the lower limit of lens redshift used in the 
calculation is 0.1, 0.2, 0.3, 0.4 and 0.5, while the upper limit is z = 1. 



of the rCDM model are the same as for SCDM except for the shape parameter (r = 0.5 for SCDM). The error in each bin is the sample 
variance corresponding to a survey area of ~ 7.84 degree 2 . The purpose of this figure is to illustrate the validity of our model for different 
cosmological models and the sensitivity of the shear kurtosis to the shape of the matter power spectrum. One can clearly see that our halo 
model reproduces the simulation results for both the shear kurtosis (upper panel) and the convergence skewness (lower panel). This success 
may be a surprise, because it has been pointed out in several works (e.g. see Figure 1 in Van Waerbeke et al. 2001b) that the JSW data does not 
exactly match the amplitude of the convergence variance predicted by the PD formula. The discrepancy may be attributed to the inaccuracy 
of the PD fitting formula for the nonlinear power spectrum. However it appears that our model can reproduce the simulation results for 
the statistical measures of weak lensing fields that are chosen to be insensitive to the power spectrum normalization but can pick up the 
non-Gaussian signals originating from the density field. Figure also shows that the skewness and kurtosis are larger for the rCDM model 
compared to the SCDM model. Thus some constraint on the shape of the power spectrum is necessary to use the non-Gaussian statistics for 
parameter estimation. 

So far our halo model calculations have assumed the mass range for the integration as 1O 3 M0 < M < 1O 18 M0. Figure [l3] plots the 
dependence of the shear kurtosis on the maximum mass cutoff used in the calculation for the SCDM model. Note that we varied the maximum 
mass cutoff for evaluations of both the shear variance and the fourth-order shear moment used in the calculation of S 7] 4; the figure shows 
the resulting dependence of S 7 ,4 on the maximum mass cutoff. It is apparent that S 7 ,4 is mainly due to massive halos with M > 1O 14 M , 
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Figure 15. In the left upper panel, the dependences of tr* (8 S ) and the shear fourth-order moment on the source galaxy redshift z s are plotted for the SCDM 
model and 9 B = 1', while the resulting dependence of the shear kurtosis is shown in the lower left panel. Note that cr* appears in the denominator of Sr/,4- 
The solid and dashed lines in the upper panel are the results for crjj and (7*), respectively, where each curve is normalized by its values for z a = 1. The right 
panel is a similar plot for the convergence skewness 5 Ki 3 given by equation ([F^l; and (k 3 ) are the denominator and numerator of S Kt :j. 



while less massive halos contribute more on smaller angular scales. More specifically, halos with M > 1O 14 M0 provide « 50% and w 99% 
contributions to the shear kurtosis at 9 S = 1' and 10', respectively. 

The dependence of the shear kurtosis on the lens redshift is shown in Figure [l4|. This figure indicates that the shear kurtosis is sensitive 
to low redshift structures with z 0.4. Next we examine the origin of the redshift dependence by plotting dependences of the numerator and 
denominator of the skewness and kurtosis separately. 

The upper-left panel in Figure [^illustrates the dependences of a® and {"ft)c on the source galaxy redshift z B for the SCDM model and 
6 S = 1' (c~ and (7*) c are the denominator and numerator of S lt 4, respectively). For comparison, the right panel shows a similar plot for 
the convergence skewness S K ,3 (the denominator and numerator of S K ,a are cr* and (k 3 ), respectively). One can see that nonlinear structures 
at lower redshifts affect the higher-order moments more strongly than the terms with powers of the variance. The lower-left panel plots the 
resulting dependence of S 7 ,4 on the source redshift and reveals that possible variation in the source redshift alters the shear kurtosis. A 
comparison of the left and right panels shows that the shear kurtosis has a stronger dependence on the source redshift than the convergence 
skewness. The shear kurtosis increases by a factor of 20 if the source redshift is varied from 0.5 to 2, while the convergence skewness varies 
about by a factor of 7. These results raise the question: what is the best survey strategy to measure the shear kurtosis? Figure [lj] suggests 
that a deeper redshift survey that probes higher-redshift structures loses some non-Gaussian signal due to projection effects. A survey to 
measure S 7 ,4 may be more efficient if it is shallower and covers greater area, provided systematic errors are well understood and the redshift 
distribution of source galaxies is known. The feasibility of the measurement of non-Gaussian statistics from lensing surveys will be presented 
in detail elsewhere (Takada et al. 2002). 

In Figure |^ we show the effect of varying the concentration parameter of the NFW profile on the shear kurtosis for the SCDM model 
and 6 S = l'. The dependence is illustrated by parameterizing the concentration parameter in terms of its normalization at the nonlinear mass 
scale today and the slope of the mass dependence as c(M, z) — co(l + z)^ 1 (M/M„ ( z — 0))~ . Here w e have again assumed that the 
redshift dependence is the same as in equation (^) as suggested by the N-body simulations (Bullock et al. 2001 ). With fixed /3, a 50% increase 
or decrease of Co leads to ~ 30% increase or decrease of 5 T ,4. Thus our results would not be strongly affected by varying c to the extent 
indicated by N-body simulations, which give a dispersion of 0.2 in lnc (Jing 2000; Bullock et al. 2001; see also Cooray & Hu 2001b for 
lensing study). On the other hand, the curves with fixed Co and varying f3 reveal that the shallower slope /3 = 0.13 leads to a larger value of 
5* 7 ,4. These results can be explained as follows. The increase of cq for a given /3 or the decrease of /3 for a given Co leads to more concentrated 
density profiles for halos more massive than the nonlinear mass scale A/*. Since these massive halos dominate the contribution to the shear 
kurtosis, this has the effect of increasing the kurtosis on the angular scales considered here. An important caveat is that the variations in 
the concentration parameter simultaneously alter the predictions for the shear variance, c?. The lower panel shows the relative errors of the 
halo model prediction to the PD results for ct 7 . It is clear that our choice (H) for c(M, z) (bold solid line) gives the closest value to the PD 
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Figure 16. In the upper panel (a), the dependence of the shear kurtosis on variations in the concentration parameter is plotted as a function of the smooth- 
ing scale. We consider the concentration parameter expressed in terms of the normalization cq at the present-day nonlinear mass scale, M*(z = 0), 
and the slope of the mass dependence as c(M, z) = cq(1 + z) {M/M t (z = 0)) - ^. The three solid lines demonstrate the dependences for (cq, 0) = 
(5, 0.2), (10, 0.2), (15, 0.2) from bottom to top, respectively, while the dashed lines are (c ,/3) = (5, 0.13), (10, 0.13), (15, 0.13). The bold solid line de- 
notes the result for (co, /3) = (10, 0.2), which we have used in this paper. In the lower panel (b), we show the relative differences between the halo model pre- 
dictions and the PD results for the shear rms. As in the upper panel, the solid lines from bottom to top are the results for (cq, /9) = (5, 0.2), (10, 0.2), (15, 0.2). 



result. Furthermore, as we have shown, our model can reproduce the simulation results for the higher-order moments of weak lensing fields. 
In this regard, therefore, as a prescription for using the halo approach to study the higher-order statistics of weak lensing, it is reasonable to 
choose the concentration parameter so that it reproduces the PD result for the variance. Conversely, if we employ a different halo profile from 
NFW, it will probably be necessary to modify our choice (||) for the concentration parameter in order to accurately describe the higher-order 
statistics. 

Finally, in Figure |lj we show the contour plot for 5* 7 ,4 with 8 S = 1' in the plane of f2 m n and as parameters for flat CDM models with 
h = 0.7. The number assigned to each contour denotes the value of log 10 5* 7 ,4, in the Sl m o — <j$ plane. Clearly, the shear kurtosis is very 
sensitive to S7 m o and has a weak dependence on ag. For example, the model with Qmo = 0.1 and as = 0.6 yields Sy,4 = 1.72 x 10 , while 
the model with f2 m o = 1.0 and the same ag leads to S" 7 , 4 = 1.46 x 10 3 . Fi gure hq shows slices of the contour plot with as = 0.6 and 
0.9, and reveals that the dependence of 5*7,4 on O m o is very close to 5* 7 ,4 oc f2~o- On the other hand, the two solid lines in Figure [l7| show 
the dependence <78f2mo\ which represents typical constraints obtained from measurements of the two-point statistics of the shear field (each 
curve is arbitrarily normalized). One can see that these two curves have a very different shape from the contours of S 7 ,4; thus measurements 
of 5 7 ,4 can break the degeneracy in the determination of f2 m o and as- Furthermore, they can cons train the dark energy c omponent of the 
universe if they are combined with the evidence for a flat universe from recent CMB measurements (Netterfield et al. 2001). 



7 DISCUSSION AND CONCLUSION 

In this paper we have investigated the kurtosis parameter of the cosmic shear field, S 7 ,4(= { r yt) c / (li) 3 ), based on the dark matter halo 
approach. The two main results revealed in this paper are summarized as follows. First, we have developed a useful approximation for 
calculating the shear kurtosis, which significantly reduces the computational time and yet provides the shear kurtosis expected within ~ 10% 
accuracy over the angular scales 1' < 9 S < 10'. Our model predictions can well match the ray-tracing simulation results for the shear kurtosis 
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Figure 17. Log contour map of the shear kurtosis parameter S 7l 4 in the Qmo-crg plane for flat CDM models with h = 0.7 and 8 B = 1'. Note that the number 
assigned to each contour is the index value of a which parameterizes the kurtosis as Sj = 10 a ; each contour is stepped by Aa = 0.2. The two solid curves 
show the dependence ctsC^q that represents typical constraints on the f2 m o-<rg plane so far obtained from the two-point shear statistics measurements (the 
two curves correspond to two different normalizations). It is clear that the curves for the two-point and four-point statistics have different shapes, and would 
therefore allow for independent determinations of the two parameters. 




Figure 18. This figure shows slices of the contour plot in Figure |lj with erg = 0.6 and 0.9, which explicitly illustrates the dependence of the shear kurtosis 
on f2 m o- The dashed line shows the dependence S lt 4 <x f2~ . 



as well as for the convergence skewness and kurtosis parameters for the SCDM, ACDM and rCDM models (see Figure pi pj, |l l| anql2|). For 
the ACDM model, the simulation data lie slightly below our predictions. It appears that the numerical results on small scales, especially for 
the higher-order moments, have not converged at the few percent level of accuracy - this is a subject that merits further investigation. While 
we have focused on lensing statistics in this paper, our results for the higher-order moments apply to the 3-dimensional density field as well. 
We show in Figure 6 that our approximations allow for the 3-dimensional skewness to be accurately computed down to sub-Mpc scales, 
which is an improvement on existing approaches in the literature. 

Second, we have shown that S 7i 4 has a strong dependence on the matter density parameter of the universe, fi m o, while it is only weakly 
dependent on the power spectrum normalization, erg, as illustrated in Figure and |T^. Thanks to this property, a measurement of Sy,4,, 
in combination with the shear two-point statistics already measured, would be valuable in constraining both fl m o and the matter power 
spectrum. For example, a marginal detection of the shear kurtosis with 50% uncertainties would yield the constraint 0.24 ,$ f2 m o £ 0.43 if 
the current concordance model with f2 m rj = 0.3, Q\o = 0.7, h = 0.7 and as — 0.9 is taken as the fiducial model. Even a null detection of 
S 7 ,4 allows us to set a lower limit on f2 m rj from the strong dependence of S 7 ,4 on low Q m o values. Thus measurements of S 7 ,4 can break the 



© 0000 RAS, MNRAS 000,[I|42^ 



The Kurtosis of the Cosmic Shear Field 2 1 



10 z 



S5 

\ 
CO 



10' 



10" 



-I 1 I 1 I 

- 


i | ' | i | i | ' | i | i |. 



AC DM 


- / \ 

N 


—2 


S 




^ \ X 










\ 








"S 

\_ \ 














: <7i 4 >o 


25 deg 2 ~~ " - . ^ ! 


■ <7 t 4 > 


100 deg 2 ^^^^ '^-^ 


i.i.i 


25 deg 2 ^^^^ 
1.1,1,1.1,1,1 



4 5 6 7 
9 [arcmin] 



10 



Figure 19. Estimates of the signal-to-noise ratio in the measurement of the connected fourth-order moment of the shear field, (jf)c, for the ACDM model 
against the smoothing scale. We assume cr e = 0.4 and n ga i = 30 arcmin -2 for the rms intrinsic ellipticity and the number density of source galaxies, and 
take the survey area to be f2 sur vcy = 25 degree 2 . For comparison, the dotted and dot-dashed lines are the results for the shear variance and ('yf)c with 



100 degree , respectively. 




degeneracies in the f2 m o and erg determination so far provided from the shear two-point statistics measurements without invoking any other 
methods. It can determine the dark energy component of the universe if combined with the strong evidence of a flat universe from the CMB 
data. It could also help resolve the puzzling inconsistency in the determination of as in the 'old' (e.g., Eke et al. 1996) and 'new' (e.g., Seljak 
2001) cluster abundance estimations (see also Van Waerbeke et al. 2002 and Lahav et al. 2001 for comments on this issue from analyses of 
the weak lensing and the galaxy redshift survey, respectively). 

We believe that the shear kurtosis is more directly applicable to data from weak lensing surveys than the well-studied higher-order 
statistics of the convergence field. To examine this issue in detail, we must examine the signal-to-noise properties of different measures of 
non-Gaussianity from realistic survey data. This would facilitate a comparison of different approaches, such as the shear kurtosis discussed 
here and the shear 3-point function proposed by Bernardeau et al (2002a). 
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For this purpose, it is crucial to correctly model possible errors in measurements of the shear, since higher-order statistics can be very 
sensitive to the noise. The main sources of error are the shot noise due to the intrinsic ellipticities of source galaxies and the sampling error 
for a finite survey area. If the intrinsic ellipticity distribution is regarded as Gaussian owing to random intrinsic orientations, we can estimate 
the dispersion for a measurement of the connected fourth-order moment of the shear field, {7*) e > following the method developed in §5 in 
Schneider et al. (1998): 



24 loV ' 



^ * T^V 5451 x + \JJ ' (43) 

where <r e is the dispersion of the intrinsic ellipticity distribution, Sl survcy the survey area and n ga i the number density of source galaxies. The 
first term on the r.h.s of the equation above denotes the sample variance and the second term the noise due to the finite number of randomly 
located source galaxy images. Here we have assumed n ga i# 2 ^> 1 and that the connected part of the higher-order moments of the shear field 
is equal to its unconnected part for simplicity pj. Figure [l^ shows an estimate of the signal-to-noise ratio in the measurement of (jf} c for 
the ACDM model at z s = 1, where we have assumed n ga i = 30 arcmin -2 and <r e = 0.4 and considered two cases of ^survey = 25 and 
100 deg 2 . Note that a signal-to-noise ratio for the measurement of S 7 ,4 is similar for the result in this figure, since the error arises mainly 
from the measurement of {"it) c compared with that of (7 2 ), as shown by comparing the solid and dotted lines. The noise is mainly due 
to the sample variance at 9 S ^ 2' for the ACDM model, while the intrinsic ellipticity noise is important at 9 S ,$,2'. One can see that for a 
survey area of 25 square degrees the measurement of {jf} c would indeed be marginally feasible on small angular scales 6 B £ 3', provided 
systematic errors can be kept under control. The results also imply an interesting possibility as discussed in Figure n5f a shallower survey, 
for a given amount of observing time, could improve the signal-to-noise ratio because the amplitude of (jt} c does not decrease as much 
as that of the variance for low redshift structures. Further the redshift distribution is easier to measure accurately for a shallower survey. In 
any case, with a survey area exceeding 100 square degrees, expected from forthcoming lensing surveys, the kurtosis measurement should be 
made with high statistical significance, as shown by the dot-dashed line, and thus prove useful for parameter estimation. 

There are some uncertainties we have ignored in the rough signal-to-noise estimate of equation (p3|). First, non-Gaussian errors are 



more important on smaller angular scales (Cooray & Hu 2001b), so sample variance must be estimated by using an adequate number of 
realizations of ray-tracing simulations or possibly by an analytical treatment using the halo approach for calculating the connected sixth- and 
eighth-order moments. Second, in actual data the noise distribution of the intrinsic ellipticities is likely to be non-Gaussian. In addition, we 
should bear in mind that the shear kurtosis is strongly affected by rare events in the tail of the measured shear distribution. Hence, we will 
need to consider some strategy to efficiently extract the shear kurtosis from realistic data which is less sensitive to unphysical rare events. 
One possible way to reduce the sample variance from such rare events is to use the probability distribution function (PDF) of the shear field, 
which is analogous to the method proposed by JSW for the study of the convergence skewness parameter. If the primordial fluctuations are 
Gaussian, the nonlinear gravitational evolution of structure formation induces the non-Gaussianity in the weak lensing fields, as investigated 
in this paper, and the weakly non-Gaussian PDF of the shear field can be modeled by the Edgeworth expansion (see Juszkiewicz et al. 1995). 
At the lowest order, we have 



P(70 = 72^) eXP (-2^ 



1 2 



(44) 



where Ha(x) = x 4 — 6x 2 + 3 is the fourth-order Hermite polynomial. The resulting PDFs for the SCDM and ACDM models are shown in 
Figure Moreover, in practice we must account for the fact that the measured shear field is a sum of the cosmic shear and noise fields. To 
obtain the PDF for the measured shear field 7° bs , therefore, we have to convolve P(7;) with the PDF of the noise, where the noise field is 
defined by smoothing the intrinsic ellipticities of source galaxies contained within top-hat apertures. Note that there are two noise fields (ei 
and 62) corresponding to the two shear fields (71 and 72), respectively. If the noise PDF is given by Pjv(ei), the PDF for the measured shear 
field, 7° bs , can be expressed by the convolution integral as 

P(% hs ) = fde Jd 1 iP N {e i )P{ rii )5 D { 1 '? s -~i i -e), (45) 

where Pjv(ei) is normalized as Jde,;Pjv(e.;) = 1. It is worth noting that the noise PDF Pjv(ei) can also be modeled in terms of the variance 
and higher-order moments of the noise field using the Edgeworth expansion. Furthermore, this method can utilize a great advantage - the 
noise PDF Pjv(ei) can be directly reconstructed from the observed shear field, e.g. by smoothing after the randomization of the position 
angle of each galaxy image. This procedure would wash out the coherent cosmic shear pattern within the smoothing aperture, but pick up 
the contribution from the intrinsic ellipticities of source galaxies. In this sense, the variance and higher-order moments of the noise field for 
a given smoothing scale can be directly extracted from the measured data, giving an estimator for the noise PDF. To obtain the noise for the 



convergence field is harder, since a non-local reconstruction from the smoothed intrinsic ellipticity field is required (Van Waerbeke 1999). 
For the shear, the theoretical model ( j45j ) for P(7° bs ) is given by a single parameter, S 7 ,4(0 S ). One can then fit the theoretical prediction to 
the measured PDF over an appropriate intermediate range of 7° bs , where the Edgeworth expansion is valid, to extract the shear kurtosis S 7l 4 
with reduced sensitivity to rare events. The quantitative improvement in the signal-to-noise will be the subject of a later study. 

Other uncertainties we have ignored in this paper are effects of the redshift distribution of source galaxies and the clustering of source 
galaxies. Figure [15] shows that the higher-order moments are more sensitive to nonlinear structures at lower redshifts than the variance and, 



" For the fourth-order moment, the ratio of the connected part to the unconnected part, (7 4 )c/[3(7 2 ) 2 ], is less than 0.5 on the angular scales we have 
considered. If this also holds for the sixth- and eighth-order moments, equation ( fl3| ) would give a conservative estimate of the noise. 
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as a result, the shear kurtosis can be sensitive to the source redshift distribution. Consider the conventionally used model for the redshift 
distribution of source galaxies 







) exp 









(46) 



with a = 2 and f3 = 1.2 (e.g., see van Waerbeke et al. (2002)) and the source redshift parameter zq — 0.48, so that it gives the mean source 
redshift (z s ) ~ 1.0. This distribution increases the value of Sj,4 at 8 S — 1' by ~ 20% for the SCDM and ACDM models, compared to the 
case with source redshifts fixed at z s =_L This increase is to some extent counterbalanced by the source clustering effect, because previous 



work based on perturbation theory (Bemardeau 1998) showed that sour ce clustering reduce s the values of the convergence skewness and 
kurtosis by ~ 10 — 20% (recettfl^rarn^rrnedj^vj^rv^tr^ 

Van Waerbeke & Mellier 2001)). It is reasonable to expect that this argument is also valid for the shear kurtosis, because we showed S~f : 4 is 
related to the convergence kurtosis S K ,4 through a geometrical constant factor (see discussions around equation (^)). 

Although the halo approach used here and in the literature assumes a spherically symmetric profile, in reality halos have non-spherical 
profiles and substructure as predicted in the CDM paradigm (e.g., Jing & Suto 2001). Our results showed that the halo model can fairly 
reproduce the ray- tracing simulation results for the 1 -point moments of the smoothed lensing fields. This success is encouraging, since the 
simulations include contributions from various realistic halo profiles. The agreement is partly because we focus only on statistical quantities 
and, therefore, to some extent the profile we need for the halo model calculation should be an average over possible halo profiles of a given 
mass. In addition, 1 -point moments of the smoothed fields are likely to be insensitive to halo profile fluctuations. We expect that the full 3- 
or 4-point correlation functions of the lensing fields would be more sensitive to profile fluctuations, because those functions should contain 
complete information on gravitational clustering up to the 3 or 4-point level through the configuration dependences. These issues will be 
presented elsewhere (Takada & Jain 2002). 

Finally, we comment on an alternative application of measurements of the shear kurtosis. If cosmological parameters including fi m o 
are precisely determined by other measurements, our results suggest that the higher-order moments of weak lensing fields could be used to 
constrain dark matter halo profiles. This problem is particularly interesting, since it can be a clue to understanding the nature of dark matter. 
Although the concentration and inner profile of dark matter halo are degenerate in giving two-point lensing statistics as argued in this paper 
(see also Seljak 2000), a detailed study may yield ways of combining the two-point and four-point shear statistics to break the degeneracy 
by exploiting the dependences of the shear kurtosis on the inner profile and the concentration parameter shown in Figure [l(| If this is the 
case, Figure |l3] indicates that measurements of the shear kurtosis at 8 B ^ 1' can constrain the properties of the halo profile at mass scales 
M > 1O 14 M . 
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benefitted from several discussions with R. Scoccimarro. We also thank E. Komatsu, A. Taruya and L. van Waerbeke for valuable discussions 
and D. Dolney for a careful reading of the manuscript. This work is supported in part by the Japan Society for Promotion of Science (JSPS) 
Research Fellowships, a NASA-LTSA grant, and a Keck foundation grant. 



APPENDIX A: PERTURBATION THEORY BISPECTRUM AND TRISPECTRUM 

The explicit forms of the bispectrum and trispectrum of the density field based on perturbation theory (e.g., Fry 1984) are 

B pt (fci,fc 2 ,fc3,2) = 2F 2 (k 1 ,k 2 )P L (k 1 ,z)P L (k 2 ,z) + 2 perm., (Al) 
T pt (fci, fe 2 , fc 3 , ki- z) = 4 [F 2 (ki 3 , -fei)F 2 (fe 13 , fc 2 )P £ (fci 3 , z)P L (k 1 , z)P L (k 2 , z) + \\ perm.] 

+6 [F 3 (fci,fc 2 ,fc 3 ) + P L {k 1 ,z)P L {k 2 ,z)P L {k z ,z) + 3 term] , (A2) 

where the redshift evolution of the linear power spectrum is given by P L (k, z) — D 2 (z)P L (k, z — 0) where D(z) is the growth factor. The 
kernels F n are calculated using perturbation theory (e.g., see Jain & Bertschniger 1994) and are expressed as 

5 l/l 1 \ „ , , 2 (fei ■ fe 2 ) 2 



F 2 (k 1 ,k 2 ) = _ + _ _ + _ (fei. fe. 



7 1 O 12 1 1.2 / V-J- ' y . 2 7. 2 

1 ju \ rvj 2 / 1 2 

7 fci2 ■ fei j . . J 1 k 2 2 (ki ■ k 2 ) 

18 fcf [F2(fe2 ' fe:5) + G2(fel ' k2)] + 18 "HfeM 



F 3 (k 1 ,k2,k 3 ) = — " 2 [F 2 (fc 2 ,fc 3 ) + G 2 (fc 1 ,fc 2 )] + — m , 2 , 2 ' [G 2 (k 2 ,k 3 ) + G 2 (k 1 ,k 2 )], (A3) 



with 



where we have ignored the extremely weak dependences of the functions F n and G 2 on cosmological parameters Qmo and f2> 



(A4) 



APPENDIX B: GEOMETRICAL PROPERTIES OF TOP-HAT FILTER FUNCTION 

The purpose of this Appendix is to derive properties of the integrals of products of the two-dimensional top-hat kernel given by equation 
(H). This is analogous to the approach in Appendix B of Bernardeau (1994). 
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We begin our discussion with deriving the following identity for the third-order products, since it is relevant for the calculation of the 
third-order moment of weak lensing fields (see 941): 



d xi f d x 2 



(2tt) 2 J (2tt) 



7 F(x 1 )F(x 2 )F(x 12 ) = 



x\dx\ 
2tt 



°X 2 dx 2 -2/ \ Wl, x 

— ( Xl > F ( X2 >> 



(Bl) 



where the top-hat kernel is F(x) = 2Ji(x)/x and xi 2 = \xi + x 2 \. 

The proof of equation (Bl ) is as follows. From the expansion formula of Bessel function (e.g., 8.532.1 in Gradshteyn & Ryzhik 2000) 
we can expand the kernel F(\xi + x 2 \) as 



F(|aJi + a3 2 |) = 2 



Ji(|xi 



■ x 2 \ 



\Xi + x 2 



4^(n+l) 



J n+ i(a;i) J n+ i(x 2 ) ryn sin[(n + 1)$ 2 ] 



(B2) 



where $2 is the angle between x\ and x 2 , giving £12 = 1 
( Jb l[ ) and then integrating it over the angle $2, we obtain 



xi x 2 sin $2 

2 + x 2 — 2xix 2 cos < I > 2 ) 1 ^ 2 . Inserting this equation into the left part of equation 



:idxi x 2 dx 2 
~2~tt 2~^~ 



F(xi)F(x 2 ) 



F( Xl )F(x 2 ) - 12 J3(ai)J2(x2) + 20^il^ + 



Xl 



X 2 



Xl 



x 2 



(B3) 



Using the recursion relation for Bessel functions and an integration formula (6.574.2 in Gradshteyn & Ryzhik 2000), the second term in the 
bracket on the r.h.s of equation above vanishes because 



f 

■Jo 



'xidxi Ji(xi) Js{xi) 



2?r 



Xl 



Xl 



°xidxi Ji(xi) 



271 



Xl 



.J 2 (xi) Ji{xi) 



Xl 



= 4-2" 1 



T(2)r(l) 



r(i)r(i/2) 



r(3/2)r(5/2)r(i/2) r(i/2)r(3/2)r(i/2) 
= 0. 



(B4) 



One can similarly find that the third term and higher terms in the bracket on the r.h.s of equation ( B3 ) vanish, and then equation ( |Bl| ) follows. 

Likewise, one can straightforwardly obtain the following identity on which the approximation for calculating the convergence forth- 
order moment is based: 



n 



d 2 Xi 
(2tt) 2 



F( Xl )F(x 2 )F(x 3 )F(x 123 ) 



n 



2tt 



F (xi)F (x 2 )F (X3), 



(B5) 



where xirj = \xi + x 2 + x 3 



APPENDIX C: 2-HALO AND 3-HALO CONTRIBUTIONS TO THE CONVERGENCE FOURTH-ORDER MOMENT 

In this Appendix, we write down the expressions for approximations used for calculations of the 2-halo and 3-halo contributions to the 
convergence fourth-order moment, which are discussed in §Q 

The 2-halo term of the convergence fourth-order moment receives two contributions, which represent taking three or two particles in 
the first halo: 



. 3 



< K 4 )f = Jd X W 4 ( X ,X,)dl 4 (x) Jt[^[ T 3i+ T i2] (h,l2,l 3 ,~li 2S )F(h)F(li 23 ), (CI) 

i — 1 

where Tjf and T^ 1 are given by equations ( p^ ) and (|l6|), and ki — li/dA(x)- It is clear that the contribution from T 3 i dominates that from 
T 22 , because the former arises mainly from 3-point correlations within one halo with a highly nonlinear density contrast, while the latter 
arises from 2-point correlations. We confirmed that the T 22 contribution is smaller even than the 3-halo term. For this reason, we ignore 
the T 22 contribution and use the following approximation for calculating the 2-halo term in the convergence fourth-order moment (we also 
include the approximation for F(h 23 ) from equation (^3[)): 



3 

,4/ \.-6/ n / ,,,,■ dn I Mi \ , . \l\dli l 2 dl 2 ~ 2 2/ 



<^ 4 )f « 4 ld X W\x,Xs)dT(x) JdMi— b(Mi,z) j^^y(h,Mi)y(l 2 ,Mi)y(li 2 ,Mi)F 2 (li)F 2 (l 2 ) 

dM2 ^2 6(M2 ' 2) (^) j l ^y^M 3 )F\h)P\h), (C2) 

where li 2 = (/ 2 + I 2 — hh) 1 ^ 2 and the factor of 4 comes from the permutation symmetry in the 2-halo trispectrum (see equation ([15])) for 
the fourth-order moment calculation. We have written the order of integration specifically to point out that the last two integrals, over M 2 
and £3, can be performed separately from the preceding three. Thus one needs to perform at most a 4-dimensional numerical integration to 

get<^ 4 ) 2h . 
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Similarly, we use the following equation to calculate the 3-halo term in the convergence fourth-order moment: 

<* 4 >f « 6 fd X W 4 (x,Xs)dJ 6 ( X ) J l ^-^BV\h,l 2 ,~h 2 )F 2 (h)F\l 2 ) f dMl ^-^b(M 1 ,z)y(l 1 ,M 1 ) 

X fajfc^KMMh,*) JdM,£- ) W) j l ^y\l 3 ,M 3 )F\ h ), (C3) 

where d 2 l 2 = l 2 ddl 2 d<& 2 , l\ 2 = (Zi + l 2 — 2l\l 2 cos $2) ly ^ 2 and the factor of 6 comes from equation (|l7|). We have used the approximation 
y(li2z) = y{h), which is valid to high accuracy for the 3-halo term. Again, a careful consideration of the dependences of the integrated 
functions reveals that we have to perform only a 4-dimensional numerical integration to obtain (k 4 )^. 
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